In this report, we reproduce the analyses using data from fMRI study 1 reported in Supplementary Material.

prep data

First, we load the relevant packages, define functions and plotting aesthetics, and load and tidy the data.

load packages

library(pacman)
pacman::p_load(tidyverse, purrr, fs, knitr, lmerTest, ggeffects, kableExtra, boot, devtools, install = TRUE)
devtools::install_github("hadley/emo")

define functions

source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")

# MLM results table function
table_model = function(model_data) {
  model_data %>%
  broom.mixed::tidy(conf.int = TRUE) %>%
  filter(effect == "fixed") %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  select(-group, -effect) %>%
  mutate_at(vars(-contains("term"), -contains("p")), round, 2) %>%
  mutate(term = gsub("cond", "", term),
         term = gsub("\\(Intercept\\)", "intercept", term),
         term = gsub("condother", "other", term),
         term = gsub("condself", "self", term),
         term = gsub("topichealth", "topic (health)", term),
         term = gsub("self_referential", "self-referential", term),
         term = gsub("self_relevance_z", "self-relevance", term),
         term = gsub("social_relevance_z", "social relevance", term),
         term = gsub(":", " x ", term),
         p = ifelse(p < .001, "< .001",
             ifelse(p > .999, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))),
         `b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", estimate, conf.low, conf.high)) %>%
  select(term, `b [95% CI]`, df, t, p)
}

simple_slopes = function(model, var, moderator, continuous = TRUE) {
  
  if (isTRUE(continuous)) {
    emmeans::emtrends(model, as.formula(paste("~", moderator)), var = var) %>%
      data.frame() %>%
      rename("trend" = 2) %>%
      mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", trend, asymp.LCL, asymp.UCL)) %>%
      select(!!moderator, `b [95% CI]`) %>%
      kable()  %>%
      kableExtra::kable_styling()
    
  } else {
    confint(emmeans::contrast(emmeans::emmeans(model, as.formula(paste("~", var, "|", moderator))), "revpairwise", by = moderator, adjust = "none")) %>%
      data.frame() %>%
      filter(grepl("control", contrast)) %>%
      mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, asymp.LCL, asymp.UCL)) %>%
      select(contrast, !!moderator, `b [95% CI]`) %>%
      arrange(contrast) %>%
      kable()  %>%
      kableExtra::kable_styling()
  }
}

define aesthetics

palette_condition = c("self" = "#ee9b00",
                      "control" = "#bb3e03",
                      "other" = "#005f73")

palette_sharing = c("#0a9396", "#ee9b00")
palette_roi = c("self-referential" = "#ee9b00",
               "mentalizing" = "#005f73")
palette_dv = c("self-relevance" = "#ee9b00",
               "social relevance" = "#005f73",
               "sharing" = "#56282D")
palette_topic = c("climate" = "#E6805E",
                 "health" = "#3A3357")

plot_aes = theme_minimal() +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        text = element_text(size = 16, family = "Futura Medium"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(color = "black"),
        axis.line = element_line(colour = "black"),
        axis.ticks.y = element_blank())

load and tidy data

merged_all = read.csv("../data/study1_data.csv")

ratings_z = merged_all %>%
  select(pID, event, trial, self_relevance, social_relevance) %>%
  unique() %>%
  mutate(self_relevance_z = scale(self_relevance, center = TRUE, scale = TRUE),
         social_relevance_z = scale(social_relevance, center = TRUE, scale = TRUE))

merged = merged_all %>%
  filter(outlier == "no" | is.na(outlier)) %>%
  group_by(pID, atlas) %>%
  mutate(parameter_estimate_std = parameter_estimate / sd(parameter_estimate, na.rm = TRUE)) %>%
  left_join(., ratings_z)

merged_wide = merged %>%
  filter(atlas %in% c("self-referential", "mentalizing")) %>%
  select(site, pID, trial, topic, cond, value, self_relevance, self_relevance_z, social_relevance, social_relevance_z, atlas, parameter_estimate_std) %>%
  spread(atlas, parameter_estimate_std) %>%
  rename("self_referential" = `self-referential`)

merged_wide_alt = merged %>%
  filter(atlas %in% c("pnas_self", "pnas_mentalizing_nopc")) %>%
  select(site, pID, trial, topic, cond, value, self_relevance, self_relevance_z, social_relevance, social_relevance_z, atlas, parameter_estimate_std) %>%
  spread(atlas, parameter_estimate_std) %>%
  rename("self_referential" = pnas_self,
         "mentalizing" = pnas_mentalizing_nopc) 



sensitivity analyses

Given the high correlation between the preregistered Neurosynth ROIs, we conducted sensitivity analyses using ROIs from Scholz et al. (2017) A neural model of valuation and information virality.

In order to maximize the differentiation between the self-referential and mentalizing ROIs, we removed the PCC/precuneus cluster from the mentalizing ROI as it overlapped with the self-referential ROI.

ROI correlations

Compared to the preregistered Neurosynth ROIs (r = .94, 95% CI [.94, .94]), the correlation between the alternative ROIs are substantially reduced.

merged_wide_alt %>%
  rmcorr::rmcorr(as.factor(pID), mentalizing, self_referential, data = .)
## 
## Repeated measures correlation
## 
## r
## 0.5588849
## 
## degrees of freedom
## 5928
## 
## p-value
## 0
## 
## 95% confidence interval
## 0.5411267 0.5761449

H1

Is greater activity in the ROIs associated with higher self and social relevance ratings?

self-referential ROI

✅ H1a: Greater activity in the self-referential ROI will be associated with higher self-relevance ratings

mod_h1a =  lmer(self_relevance ~ self_referential + (1 + self_referential | pID),
               data = merged_wide_alt,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h1a = table_model(mod_h1a)

table_h1a %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.56 [2.48, 2.64] 84.62 65.06 < .001
self-referential 0.03 [0.01, 0.06] 83.66 2.40 .018

summary

summary(mod_h1a)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: self_relevance ~ self_referential + (1 + self_referential | pID)
##    Data: merged_wide_alt
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 16758.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4707 -0.7086  0.1383  0.6783  2.4399 
## 
## Random effects:
##  Groups   Name             Variance Std.Dev. Corr 
##  pID      (Intercept)      0.116790 0.34175       
##           self_referential 0.001606 0.04008  -0.79
##  Residual                  0.918028 0.95814       
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##                  Estimate Std. Error       df t value            Pr(>|t|)    
## (Intercept)       2.55698    0.03930 84.61525  65.056 <0.0000000000000002 ***
## self_referential  0.03131    0.01303 83.66472   2.403              0.0185 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## self_rfrntl -0.344

mentalizing ROI

✅ H1b: Greater activity in the mentalizing ROI will be associated with higher social relevance ratings

mod_h1b = lmer(social_relevance ~ mentalizing + (1 + mentalizing | pID),
               data = merged_wide_alt,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h1b = table_model(mod_h1b)

table_h1b %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.67 [2.58, 2.75] 84.07 64.07 < .001
mentalizing 0.03 [0.01, 0.06] 82.98 2.86 .005

summary

summary(mod_h1b)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: social_relevance ~ mentalizing + (1 + mentalizing | pID)
##    Data: merged_wide_alt
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 15841.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8050 -0.7206  0.1710  0.6480  2.6691 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  pID      (Intercept) 0.134851 0.36722       
##           mentalizing 0.001039 0.03224  -0.09
##  Residual             0.784321 0.88562       
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)  2.66512    0.04160 84.07281  64.067 < 0.0000000000000002 ***
## mentalizing  0.03406    0.01192 82.97654   2.856              0.00541 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## mentalizing -0.100

combined plot

vals = seq(-4.5,4.5,.1)
predicted_h1 = ggeffects::ggpredict(mod_h1a, c("self_referential [vals]")) %>%
  data.frame() %>%
  mutate(roi = "self-referential",
         variable = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h1b, c("mentalizing [vals]")) %>%
              data.frame() %>%
              mutate(roi = "mentalizing",
                     variable = "social relevance"))

predicted_sub_h1 = ggeffects::ggpredict(mod_h1a, terms = c("self_referential [vals]", "pID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "self-referential",
         variable = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h1b, c("mentalizing [vals]", "pID"), type = "random") %>%
              data.frame() %>%
              mutate(roi = "mentalizing",
                     variable = "social relevance"))

predicted_h1 %>%
  ggplot(aes(x, predicted)) +
  stat_smooth(data = predicted_sub_h1, aes(group = group, color = roi), geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = roi), alpha = .5, color = NA) +
  geom_line(aes(color = roi), size = 2) +
  facet_grid(~variable) +
  scale_color_manual(name = "", values = palette_roi, guide = FALSE) +
  scale_fill_manual(name = "", values = palette_roi, guide = FALSE) +
  labs(x = "\nROI activity (SD)", y = "predicted rating\n") +
  plot_aes

H4

Do the manipulations increase neural activity in brain regions associated with self-referential processing and mentalizing?

self-referential ROI

✅ H4a: Self-focused intervention (compared to control) will increase brain activity in ROIs related to self-referential processes.

mod_h4a = lmer(self_referential ~ cond + (1 + cond | pID),
              data = merged_wide_alt,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h4a = table_model(mod_h4a)

table_h4a %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 0.23 [0.12, 0.33] 84.12 4.38 < .001
other 0.11 [0.03, 0.19] 84.42 2.86 .005
self 0.13 [0.04, 0.21] 83.69 2.84 .006

summary

summary(mod_h4a)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: self_referential ~ cond + (1 + cond | pID)
##    Data: merged_wide_alt
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 17254.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4678 -0.6543 -0.0102  0.6474  3.5891 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  pID      (Intercept) 0.18879  0.4345              
##           condother   0.04997  0.2235   -0.05      
##           condself    0.08492  0.2914    0.05  0.49
##  Residual             0.97421  0.9870              
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##             Estimate Std. Error       df t value  Pr(>|t|)    
## (Intercept)  0.22817    0.05204 84.12046   4.385 0.0000334 ***
## condother    0.11309    0.03950 84.41676   2.863   0.00529 ** 
## condself     0.12606    0.04441 83.69460   2.838   0.00569 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) cndthr
## condother -0.267       
## condself  -0.176  0.493

mentalizing ROI

❌ H4b: Other-focused intervention (compared to control) will increase brain activity in ROIs related to mentalizing processes.

mod_h4b = lmer(mentalizing ~ cond + (1 + cond | pID),
              data = merged_wide_alt,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h4b = table_model(mod_h4b)

table_h4b %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 0.27 [0.15, 0.39] 84.19 4.58 < .001
other 0.01 [-0.07, 0.08] 83.16 0.21 .837
self 0.05 [-0.04, 0.13] 84.03 1.11 .271

summary

summary(mod_h4b)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mentalizing ~ cond + (1 + cond | pID)
##    Data: merged_wide_alt
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 17292.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4694 -0.6574 -0.0019  0.6645  4.0424 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  pID      (Intercept) 0.25179  0.5018              
##           condother   0.04215  0.2053   -0.45      
##           condself    0.07355  0.2712   -0.25  0.87
##  Residual             0.98487  0.9924              
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  0.269101   0.058773 84.192890   4.579 0.000016 ***
## condother    0.007961   0.038457 83.160115   0.207    0.837    
## condself     0.047633   0.043001 84.025816   1.108    0.271    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) cndthr
## condother -0.461       
## condself  -0.354  0.641

combined plot

predicted_h4 = ggeffects::ggpredict(mod_h4a, c("cond")) %>%
  data.frame() %>%
  mutate(atlas = "self-referential") %>%
  bind_rows(ggeffects::ggpredict(mod_h4b, c("cond")) %>%
              data.frame() %>%
              mutate(atlas = "mentalizing")) %>%
  mutate(x = factor(x, levels = c("self", "control", "other")),
         atlas = factor(atlas, levels = c("self-referential", "mentalizing")))

predicted_sub_h4 = ggeffects::ggpredict(mod_h4a, terms = c("cond", "pID"), type = "random") %>%
  data.frame() %>%
  mutate(atlas = "self-referential") %>%
  bind_rows(ggeffects::ggpredict(mod_h4b, c("cond", "pID"), type = "random") %>%
              data.frame() %>%
              mutate(atlas = "mentalizing")) %>%
  mutate(x = factor(x, levels = c("self", "control", "other")),
         atlas = factor(atlas, levels = c("self-referential", "mentalizing")))

predicted_h4 %>%
  ggplot(aes(x = x, y = predicted)) +
  stat_summary(data = predicted_sub_h4, aes(group = group), fun = "mean", geom = "line",
               size = .1, color = "grey50") +
  stat_summary(aes(group = group), fun = "mean", geom = "line", size = 1) +
  geom_pointrange(aes(color = x, ymin = conf.low, ymax = conf.high), size = .75) +
  facet_grid(~atlas) +
  scale_color_manual(name = "", values = palette_condition, guide = "none") +
  scale_alpha_manual(name = "", values = c(1, .5)) +
  labs(x = "", y = "ROI activity (SD)\n") +
  plot_aes

H6

Is ROI activity positively related to sharing intentions?

self-referential ROI

✅ Stronger activity in the self-referential ROI will be related to higher sharing intentions.

mod_h6a = lmer(value ~ self_referential + (1 + self_referential | pID),
              data = merged_wide_alt,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h6a = table_model(mod_h6a)

table_h6a %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.60 [2.52, 2.68] 85.08 67.91 < .001
self-referential 0.06 [0.03, 0.09] 83.56 4.40 < .001

summary

summary(mod_h6a)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: value ~ self_referential + (1 + self_referential | pID)
##    Data: merged_wide_alt
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 16641.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5502 -0.7218  0.1157  0.7344  2.1951 
## 
## Random effects:
##  Groups   Name             Variance Std.Dev. Corr 
##  pID      (Intercept)      0.10916  0.33040       
##           self_referential 0.00401  0.06333  -0.36
##  Residual                  0.93217  0.96549       
## Number of obs: 5935, groups:  pID, 85
## 
## Fixed effects:
##                  Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)       2.59963    0.03828 85.07741  67.910 < 0.0000000000000002 ***
## self_referential  0.06272    0.01427 83.55508   4.396            0.0000323 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## self_rfrntl -0.251

mentalizing ROI

✅ Stronger activation in the mentalizing ROI will be related to higher sharing intentions.

mod_h6b = lmer(value ~ mentalizing + (1 + mentalizing | pID),
              data = merged_wide_alt,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h6b = table_model(mod_h6b)

table_h6b %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.61 [2.53, 2.68] 84.71 68.68 < .001
mentalizing 0.04 [0.02, 0.07] 83.49 3.41 .001

summary

summary(mod_h6b)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: value ~ mentalizing + (1 + mentalizing | pID)
##    Data: merged_wide_alt
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 16658.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6298 -0.7135  0.1138  0.7344  2.0984 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  pID      (Intercept) 0.107678 0.32814      
##           mentalizing 0.001175 0.03428  0.08
##  Residual             0.936640 0.96780      
## Number of obs: 5935, groups:  pID, 85
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)  2.60723    0.03796 84.70738  68.683 < 0.0000000000000002 ***
## mentalizing  0.04450    0.01306 83.48501   3.407              0.00101 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## mentalizing -0.070

combined plot

vals = seq(-4.5,4.5,.1)

predicted_h6 = ggeffects::ggpredict(mod_h6a, c("self_referential [vals]")) %>%
  data.frame() %>%
  mutate(roi = "self-referential") %>%
  bind_rows(ggeffects::ggpredict(mod_h6b, c("mentalizing [vals]")) %>%
              data.frame() %>%
              mutate(roi = "mentalizing")) %>%
  mutate(roi = factor(roi, levels = c("self-referential", "mentalizing")))

predicted_sub_h6 = ggeffects::ggpredict(mod_h6a, terms = c("self_referential [vals]", "pID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "self-referential") %>%
  bind_rows(ggeffects::ggpredict(mod_h6b, c("mentalizing [vals]", "pID"), type = "random") %>%
              data.frame() %>%
              mutate(roi = "mentalizing")) %>%
  mutate(roi = factor(roi, levels = c("self-referential", "mentalizing")))

predicted_h6 %>%
  ggplot(aes(x = x, y = predicted, color = roi, fill = roi)) +
  stat_smooth(data = predicted_sub_h6, aes(group = group), geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .2, color = NA) +
  geom_line(size = 2) +
  facet_grid(~roi) +
  scale_color_manual(name = "", values = palette_roi) +
  scale_fill_manual(name = "", values = palette_roi) +
  labs(y = "predicted sharing intention\n", x = "\nROI activity (SD)") +
  plot_aes +
  theme(legend.position = "none")

H7

Is there an indirect effect of the condition on sharing intentions through activity in self-referential and mentalizing ROIs?

prep data

# source functions
source("indirectMLM.R")

# create self condition dataframe
data_med_self = merged_wide_alt %>%
  filter(!cond == "other") %>%
  mutate(cond = ifelse(cond == "self", 1, 0)) %>%
  select(pID, site, trial, cond, value, self_referential) %>%
  data.frame()

# create social condition dataframe
data_med_other = merged_wide_alt %>%
  filter(!cond == "self") %>%
  mutate(cond = ifelse(cond == "other", 1, 0)) %>%
  select(pID, site, trial, cond, value, mentalizing) %>%
  data.frame()

# define variables
y_var = "value"

self condition

✅ H7a: The effect of self-focused intervention on sharing intention will be mediated by increased activity in the self-referential ROI.

model_name = "mediation_self"
data = data_med_self

if (file.exists(sprintf("models/model_%s_alternative.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("models/model_%s_alternative.RDS", model_name)))
} else {
  assign(get("model_name"), boot(data = data, statistic = indirect.mlm, R = 500,
                                 y = y_var, x = "cond", mediator = "self_referential", group.id = "pID",
                                 between.m = F, uncentered.x = F))
  saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s_alternative.RDS", model_name))
}

indirect.mlm.summary(get(model_name))
## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.006, 0.009]
## 
## 
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: 0.009 [0.001, 0.019]
## Biased Estimate of Within-subjects Indirect Effect: 0.009 [0.004, 0.015]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.009]
## 
## 
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.046 [-0.106, 0.013]
## Biased Total Effect of X on Y (c path): -0.044 [-0.105, 0.016]
## Bias in Total Effect: 0.002 [0, 0.006]
## 
## 
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.055 [-0.12, 0.002]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): 0.126 [0.064, 0.184]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.069 [0.041, 0.1]

other condition

❌ H7b: The effect of other-focused intervention on sharing intention will be mediated by increased activity in the mentalizing ROI.

model_name = "mediation_other"
data = data_med_other

if (file.exists(sprintf("models/model_%s_alternative.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("models/model_%s_alternative.RDS", model_name)))
} else {
  assign(get("model_name"), boot(data = data, statistic = indirect.mlm, R = 500,
                                 y = y_var, x = "cond", mediator = "mentalizing", group.id = "pID",
                                 between.m = F, uncentered.x = F))
  saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s_alternative.RDS", model_name))
}

indirect.mlm.summary(get(model_name))
## #### Population Covariance ####
## Covariance of Random Slopes a and b: -0.001 [-0.005, 0.006]
## 
## 
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: 0 [-0.005, 0.007]
## Biased Estimate of Within-subjects Indirect Effect: 0 [-0.003, 0.004]
## Bias in Within-subjects Indirect Effect: 0.001 [0, 0.006]
## 
## 
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.031 [-0.089, 0.028]
## Biased Total Effect of X on Y (c path): -0.032 [-0.09, 0.028]
## Bias in Total Effect: 0.002 [0, 0.005]
## 
## 
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.03 [-0.091, 0.027]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): 0.008 [-0.053, 0.067]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.052 [0.023, 0.084]



combined table

table_h1a %>% mutate(DV = "H1a: Self-relevance") %>%
  bind_rows(table_h1b %>% mutate(DV = "H1b: Social relevance")) %>%
  bind_rows(table_h4a %>% mutate(DV = "H4a: Self-referential ROI")) %>%
  bind_rows(table_h4b %>% mutate(DV = "H4b: Mentalizing ROI")) %>%
  bind_rows(table_h6a %>% mutate(DV = "H6a: Sharing intention")) %>%
  bind_rows(table_h6b %>% mutate(DV = "H6b: Sharing intention")) %>%
  select(DV, everything()) %>%
  kable() %>%
  kable_styling()
DV term b [95% CI] df t p
H1a: Self-relevance intercept 2.56 [2.48, 2.64] 84.62 65.06 < .001
H1a: Self-relevance self-referential 0.03 [0.01, 0.06] 83.66 2.40 .018
H1b: Social relevance intercept 2.67 [2.58, 2.75] 84.07 64.07 < .001
H1b: Social relevance mentalizing 0.03 [0.01, 0.06] 82.98 2.86 .005
H4a: Self-referential ROI intercept 0.23 [0.12, 0.33] 84.12 4.38 < .001
H4a: Self-referential ROI other 0.11 [0.03, 0.19] 84.42 2.86 .005
H4a: Self-referential ROI self 0.13 [0.04, 0.21] 83.69 2.84 .006
H4b: Mentalizing ROI intercept 0.27 [0.15, 0.39] 84.19 4.58 < .001
H4b: Mentalizing ROI other 0.01 [-0.07, 0.08] 83.16 0.21 .837
H4b: Mentalizing ROI self 0.05 [-0.04, 0.13] 84.03 1.11 .271
H6a: Sharing intention intercept 2.60 [2.52, 2.68] 85.08 67.91 < .001
H6a: Sharing intention self-referential 0.06 [0.03, 0.09] 83.56 4.40 < .001
H6b: Sharing intention intercept 2.61 [2.53, 2.68] 84.71 68.68 < .001
H6b: Sharing intention mentalizing 0.04 [0.02, 0.07] 83.49 3.41 .001

moderation by article topic

These analyses explore whether the analyses reported in study 1 of the main manuscript are moderated by article topic (health or climate).

H1

Are the relationships between ROI activity and self and social relevance ratings moderated by article topic?

self-referential ROI

There is a main effect of topic, such that health articles elicited greater activity in the self-referential ROI compared to climate articles.

These data are not consistent with moderation by topic.

mod_h1am =  lmer(self_relevance ~ self_referential * topic + (1 + topic | pID),
               data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h1am = table_model(mod_h1am)

table_h1am %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.49 [2.40, 2.59] 83.88 51.73 < .001
self-referential 0.03 [-0.01, 0.06] 5807.44 1.56 .120
topic (health) 0.13 [0.01, 0.26] 84.45 2.14 .035
self-referential x topic (health) 0.04 [-0.01, 0.08] 5759.32 1.49 .137

summary

summary(mod_h1am)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: self_relevance ~ self_referential * topic + (1 + topic | pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 16417.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.93065 -0.70980  0.09214  0.69503  2.84754 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  pID      (Intercept) 0.1735   0.4166        
##           topichealth 0.2838   0.5327   -0.59
##  Residual             0.8424   0.9178        
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##                                Estimate Std. Error         df t value
## (Intercept)                     2.49412    0.04821   83.87912  51.731
## self_referential                0.02609    0.01676 5807.44003   1.557
## topichealth                     0.13385    0.06255   84.45350   2.140
## self_referential:topichealth    0.03518    0.02364 5759.31664   1.488
##                                         Pr(>|t|)    
## (Intercept)                  <0.0000000000000002 ***
## self_referential                          0.1195    
## topichealth                               0.0353 *  
## self_referential:topichealth              0.1368    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) slf_rf tpchlt
## self_rfrntl -0.030              
## topichealth -0.606  0.023       
## slf_rfrntl:  0.021 -0.705 -0.053

mentalizing ROI

There is a main effect of topic, such that health articles elicited greater activity in the mentalizing ROI compared to climate articles.

These data are not consistent with moderation by topic.

mod_h1bm = lmer(social_relevance ~ mentalizing * topic + (1 + mentalizing | pID),
               data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h1bm = table_model(mod_h1bm)

table_h1bm %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.51 [2.43, 2.60] 97.93 58.12 < .001
mentalizing 0.03 [0.00, 0.07] 226.60 2.16 .032
topic (health) 0.29 [0.24, 0.34] 5924.99 12.18 < .001
mentalizing x topic (health) 0.01 [-0.03, 0.05] 5903.41 0.66 .512

summary

summary(mod_h1bm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: social_relevance ~ mentalizing * topic + (1 + mentalizing | pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 15673.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7740 -0.7047  0.1262  0.6762  2.7725 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  pID      (Intercept) 0.134996 0.36742       
##           mentalizing 0.001968 0.04436  -0.17
##  Residual             0.760437 0.87203       
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##                           Estimate Std. Error         df t value
## (Intercept)                2.51344    0.04325   97.92856  58.119
## mentalizing                0.03426    0.01585  226.60228   2.162
## topichealth                0.29013    0.02382 5924.99396  12.178
## mentalizing:topichealth    0.01347    0.02055 5903.41026   0.655
##                                    Pr(>|t|)    
## (Intercept)             <0.0000000000000002 ***
## mentalizing                          0.0317 *  
## topichealth             <0.0000000000000002 ***
## mentalizing:topichealth              0.5123    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mntlzn tpchlt
## mentalizing -0.156              
## topichealth -0.267  0.173       
## mntlzng:tpc  0.076 -0.636 -0.323

combined plot

vals = seq(-4.5,4.5,.1)

predicted_h1m = ggeffects::ggpredict(mod_h1am, c("self_referential [vals]", "topic")) %>%
  data.frame() %>%
  mutate(roi = "self-referential",
         variable = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h1bm, c("mentalizing [vals]", "topic")) %>%
              data.frame() %>%
              mutate(roi = "mentalizing",
                     variable = "social relevance"))


predicted_sub_h1m = ggeffects::ggpredict(mod_h1am, terms = c("self_referential [vals]", "topic", "pID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "self-referential",
         variable = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h1bm, c("mentalizing [vals]", "topic", "pID"), type = "random") %>%
              data.frame() %>%
              mutate(roi = "mentalizing",
                     variable = "social relevance")) 

predicted_h1m %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  stat_smooth(data = predicted_sub_h1m, aes(group = interaction(group, facet)),
              geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, color = NA) +
  geom_line(size = 2) +
  facet_grid(~variable) +
  scale_color_manual(name = "", values = palette_topic) +
  scale_fill_manual(name = "", values = palette_topic) +
  labs(x = "\nROI activity (SD)", y = "predicted rating\n") +
  plot_aes +
  theme(legend.key.width=unit(2,"cm"))

H2

Are the effects of the experimental manipulations on relevance moderated by article topic?

self-relevance

There is a main effect of topic such that health articles are rated as more self-relevant than climate articles.

The was also an interaction such that the effect of the self-focused condition on self-relevance was weaker for health articles.

mod_h2am = lmer(self_relevance ~ cond * topic + (1 + topic | pID),
               data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h2am = table_model(mod_h2am)

table_h2am %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.44 [2.34, 2.55] 128.84 45.55 < .001
other 0.04 [-0.04, 0.12] 5840.55 1.03 .304
self 0.12 [0.04, 0.20] 5840.64 2.84 .004
topic (health) 0.22 [0.08, 0.36] 139.35 3.18 .002
other x topic (health) -0.07 [-0.19, 0.04] 5840.57 -1.28 .202
self x topic (health) -0.17 [-0.28, -0.06] 5840.67 -2.92 .004

summary

summary(mod_h2am)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: self_relevance ~ cond * topic + (1 + topic | pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 16429.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.86770 -0.70174  0.07785  0.70050  2.75230 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  pID      (Intercept) 0.1732   0.4162        
##           topichealth 0.2825   0.5315   -0.59
##  Residual             0.8436   0.9185        
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##                         Estimate Std. Error         df t value
## (Intercept)              2.44347    0.05364  128.84030  45.551
## condother                0.04219    0.04103 5840.55170   1.028
## condself                 0.11670    0.04102 5840.64260   2.845
## topichealth              0.22474    0.07077  139.35146   3.176
## condother:topichealth   -0.07404    0.05803 5840.57270  -1.276
## condself:topichealth    -0.16935    0.05805 5840.67318  -2.918
##                                   Pr(>|t|)    
## (Intercept)           < 0.0000000000000002 ***
## condother                          0.30388    
## condself                           0.00446 ** 
## topichealth                        0.00184 ** 
## condother:topichealth              0.20208    
## condself:topichealth               0.00354 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndthr cndslf tpchlt cndth:
## condother   -0.381                            
## condself    -0.382  0.499                     
## topichealth -0.628  0.289  0.289              
## cndthr:tpch  0.270 -0.707 -0.353 -0.410       
## cndslf:tpch  0.270 -0.353 -0.707 -0.410  0.500

social relevance

There is a main effect of topic such that health articles are rated as more socially relevant than climate articles.

These data are not consistent with moderation by topic.

mod_h2bm = lmer(social_relevance ~ cond * topic + (1 + topic | pID),
               data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h2bm = table_model(mod_h2bm)

table_h2bm %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.49 [2.39, 2.60] 122.27 47.37 < .001
other 0.03 [-0.04, 0.10] 5840.82 0.78 .437
self 0.07 [-0.00, 0.15] 5840.90 1.95 .051
topic (health) 0.31 [0.19, 0.43] 151.89 5.04 < .001
other x topic (health) 0.03 [-0.07, 0.14] 5840.82 0.58 .561
self x topic (health) -0.06 [-0.16, 0.05] 5840.94 -1.08 .282

summary

summary(mod_h2bm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: social_relevance ~ cond * topic + (1 + topic | pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 15446.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9667 -0.6908  0.1057  0.6750  3.0215 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  pID      (Intercept) 0.1747   0.4179        
##           topichealth 0.1938   0.4402   -0.49
##  Residual             0.7150   0.8456        
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##                         Estimate Std. Error         df t value
## (Intercept)              2.49147    0.05260  122.26926  47.368
## condother                0.02939    0.03778 5840.81532   0.778
## condself                 0.07359    0.03777 5840.90467   1.949
## topichealth              0.30674    0.06089  151.88841   5.037
## condother:topichealth    0.03107    0.05342 5840.82084   0.582
## condself:topichealth    -0.05745    0.05344 5840.93728  -1.075
##                                   Pr(>|t|)    
## (Intercept)           < 0.0000000000000002 ***
## condother                           0.4366    
## condself                            0.0514 .  
## topichealth                     0.00000132 ***
## condother:topichealth               0.5608    
## condself:topichealth                0.2824    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndthr cndslf tpchlt cndth:
## condother   -0.358                            
## condself    -0.358  0.499                     
## topichealth -0.552  0.309  0.309              
## cndthr:tpch  0.253 -0.707 -0.353 -0.439       
## cndslf:tpch  0.253 -0.353 -0.707 -0.439  0.500

combined plot

predicted_h2m = ggeffects::ggpredict(mod_h2am, c("cond", "topic")) %>%
  data.frame() %>%
  mutate(model = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h2bm, c("cond", "topic")) %>%
              data.frame() %>%
              mutate(model = "social relevance")) %>%
  mutate(x = factor(x, levels = c("self", "control", "other")))

predicted_sub_h2m = ggeffects::ggpredict(mod_h2am, terms = c("cond", "topic", "pID"), type = "random") %>%
  data.frame() %>%
  mutate(model = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h2bm, c("cond", "topic", "pID"), type = "random") %>%
              data.frame() %>%
              mutate(model = "social relevance")) %>%
  mutate(x = factor(x, levels = c("self", "control", "other"))) 
  
predicted_h2m %>%
  ggplot(aes(x = x, y = predicted, color = group)) +
  stat_summary(data = predicted_sub_h2m, aes(group = interaction(group, facet)), fun = "mean", geom = "line", size = .1, alpha = .5) +
  stat_summary(aes(group = group), fun = "mean", geom = "line", size = 1) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high, group = group),
                  size = .75) +
  facet_grid(~model) +
  scale_color_manual(name = "", values = palette_topic) +
  labs(x = "", y = "predicted rating\n") +
  plot_aes +
  theme(legend.position = c(.85, .15))

H3

Are the relationships between self and social relevance and sharing intentions moderated by article topic?

The relationship between self-relevance and sharing intentions was not moderated by topic.

However, the relationship between social relevance and sharing intentions was slightly stronger for health articles compared to climate articles.

mod_h3m = lmer(value ~ self_relevance_z * topic + social_relevance_z * topic + (1 + self_relevance_z + social_relevance_z | pID),
              data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

plot

predicted_h3m = ggeffects::ggpredict(mod_h3m, c("self_relevance_z", "topic")) %>%
  data.frame() %>%
  mutate(variable = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h3m, c("social_relevance_z", "topic")) %>%
              data.frame() %>%
              mutate(variable = "social relevance"))

predicted_sub_h3m = ggeffects::ggpredict(mod_h3m, terms = c("self_relevance_z", "topic", "pID"), type = "random") %>%
  data.frame() %>%
  mutate(variable = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h3m, c("social_relevance_z", "topic", "pID"), type = "random") %>%
              data.frame() %>%
              mutate(variable = "social relevance"))

predicted_h3m %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  stat_smooth(data = predicted_sub_h3m, aes(group = interaction(group, facet)),
              geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, color = NA) +
  geom_line(size = 2) +
  facet_grid(~variable) +
  scale_color_manual(name = "", values = palette_topic) +
  scale_fill_manual(name = "", values = palette_topic) +
  labs(x = "\nrating (SD)", y = "predicted sharing intention\n") +
  plot_aes +
  theme(legend.key.width=unit(2,"cm"))

model table

table_h3m = table_model(mod_h3m)

table_h3m %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.56 [2.50, 2.61] 113.93 86.74 < .001
self-relevance 0.30 [0.24, 0.35] 275.19 10.76 < .001
topic (health) 0.14 [0.10, 0.19] 5544.95 6.43 < .001
social relevance 0.20 [0.14, 0.26] 180.07 6.53 < .001
self-relevance x topic (health) 0.03 [-0.03, 0.09] 5555.76 1.04 .297
topic (health) x social relevance 0.06 [0.00, 0.12] 5628.34 1.99 .047

summary

summary(mod_h3m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: value ~ self_relevance_z * topic + social_relevance_z * topic +  
##     (1 + self_relevance_z + social_relevance_z | pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 14865.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2186 -0.6964  0.0530  0.6918  3.0662 
## 
## Random effects:
##  Groups   Name               Variance Std.Dev. Corr       
##  pID      (Intercept)        0.05174  0.2275              
##           self_relevance_z   0.01335  0.1155   -0.40      
##           social_relevance_z 0.02908  0.1705    0.20 -0.52
##  Residual                    0.68009  0.8247              
## Number of obs: 5935, groups:  pID, 85
## 
## Fixed effects:
##                                  Estimate Std. Error         df t value
## (Intercept)                       2.55602    0.02947  113.92813  86.741
## self_relevance_z                  0.29550    0.02746  275.19097  10.761
## topichealth                       0.14413    0.02243 5544.94561   6.426
## social_relevance_z                0.19768    0.03027  180.06605   6.530
## self_relevance_z:topichealth      0.03216    0.03081 5555.75850   1.044
## topichealth:social_relevance_z    0.06087    0.03064 5628.34238   1.986
##                                            Pr(>|t|)    
## (Intercept)                    < 0.0000000000000002 ***
## self_relevance_z               < 0.0000000000000002 ***
## topichealth                          0.000000000142 ***
## social_relevance_z                   0.000000000648 ***
## self_relevance_z:topichealth                  0.297    
## topichealth:social_relevance_z                0.047 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) slf_r_ tpchlt scl_r_ slf__:
## slf_rlvnc_z -0.180                            
## topichealth -0.378  0.032                     
## scl_rlvnc_z  0.160 -0.655 -0.081              
## slf_rlvnc_:  0.019 -0.679 -0.020  0.435       
## tpchlth:s__ -0.053  0.477  0.006 -0.575 -0.675

H4

Are the effects of the experimental manipulations on ROI activity moderated by article topic?

self-referential ROI

There is a main effect of topic, such that health articles elicited greater activity in the self-referential ROI compared to climate articles.

These data are not consistent with moderation by topic.

mod_h4am = lmer(self_referential ~ cond * topic + (1 + cond + topic | pID),
              data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h4am = table_model(mod_h4am)

table_h4am %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 0.01 [-0.10, 0.13] 101.71 0.24 .812
other 0.09 [-0.01, 0.19] 222.78 1.76 .080
self 0.13 [0.02, 0.23] 193.75 2.38 .018
topic (health) 0.14 [0.05, 0.23] 528.65 2.99 .003
other x topic (health) -0.01 [-0.13, 0.12] 5677.79 -0.09 .928
self x topic (health) -0.08 [-0.20, 0.05] 5678.30 -1.22 .221

summary

summary(mod_h4am)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: self_referential ~ cond * topic + (1 + cond + topic | pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 17272.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8707 -0.6542  0.0020  0.6446  3.5931 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr             
##  pID      (Intercept) 0.21648  0.4653                    
##           condother   0.04666  0.2160   -0.16            
##           condself    0.07433  0.2726   -0.06  0.59      
##           topichealth 0.01628  0.1276    0.21 -0.25 -0.10
##  Residual             0.97242  0.9861                    
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##                          Estimate  Std. Error          df t value Pr(>|t|)   
## (Intercept)              0.014128    0.059286  101.710846   0.238  0.81212   
## condother                0.087855    0.049901  222.781895   1.761  0.07968 . 
## condself                 0.126311    0.053055  193.748663   2.381  0.01825 * 
## topichealth              0.138255    0.046191  528.653661   2.993  0.00289 **
## condother:topichealth   -0.005645    0.062304 5677.787287  -0.091  0.92781   
## condself:topichealth    -0.076332    0.062321 5678.302007  -1.225  0.22069   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndthr cndslf tpchlt cndth:
## condother   -0.390                            
## condself    -0.337  0.520                     
## topichealth -0.299  0.385  0.378              
## cndthr:tpch  0.262 -0.624 -0.293 -0.675       
## cndslf:tpch  0.262 -0.311 -0.587 -0.675  0.500

mentalizing ROI

There is a main effect of topic, such that health articles elicited greater activity in the mentalizing ROI compared to climate articles.

These data are not consistent with moderation by topic.

mod_h4bm = lmer(mentalizing ~ cond * topic + (1 + cond + topic | pID),
              data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h4bm = table_model(mod_h4bm)

table_h4bm %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 0.27 [0.16, 0.39] 102.44 4.63 < .001
other 0.06 [-0.04, 0.15] 233.15 1.16 .249
self 0.10 [-0.00, 0.20] 198.21 1.91 .057
topic (health) 0.11 [0.03, 0.20] 606.69 2.52 .012
other x topic (health) 0.01 [-0.12, 0.13] 5677.90 0.09 .927
self x topic (health) -0.06 [-0.18, 0.07] 5678.40 -0.89 .373

summary

summary(mod_h4bm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mentalizing ~ cond * topic + (1 + cond + topic | pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 17283.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6798 -0.6607  0.0161  0.6711  3.2619 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr             
##  pID      (Intercept) 0.208790 0.45694                   
##           condother   0.039235 0.19808  -0.14            
##           condself    0.069761 0.26412  -0.04  0.61      
##           topichealth 0.008325 0.09124   0.19 -0.51 -0.17
##  Residual             0.977955 0.98892                   
## Number of obs: 6014, groups:  pID, 85
## 
## Fixed effects:
##                         Estimate Std. Error         df t value  Pr(>|t|)    
## (Intercept)              0.27121    0.05856  102.43949   4.631 0.0000107 ***
## condother                0.05683    0.04913  233.15144   1.157    0.2486    
## condself                 0.10063    0.05265  198.21060   1.911    0.0574 .  
## topichealth              0.11419    0.04529  606.68900   2.521    0.0119 *  
## condother:topichealth    0.00576    0.06248 5677.89718   0.092    0.9265    
## condself:topichealth    -0.05564    0.06250 5678.40382  -0.890    0.3734    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndthr cndslf tpchlt cndth:
## condother   -0.391                            
## condself    -0.332  0.521                     
## topichealth -0.332  0.389  0.388              
## cndthr:tpch  0.266 -0.636 -0.296 -0.690       
## cndslf:tpch  0.266 -0.317 -0.593 -0.690  0.500

combined plot

predicted_h4m = ggeffects::ggpredict(mod_h4am, c("cond", "topic")) %>%
  data.frame() %>%
  mutate(atlas = "self-referential") %>%
  bind_rows(ggeffects::ggpredict(mod_h4bm, c("cond", "topic")) %>%
              data.frame() %>%
              mutate(atlas = "mentalizing")) %>%
  mutate(x = factor(x, levels = c("self", "control", "other")),
         atlas = factor(atlas, levels = c("self-referential", "mentalizing")))

predicted_sub_h4m = ggeffects::ggpredict(mod_h4am, terms = c("cond", "topic", "pID"), type = "random") %>%
  data.frame() %>%
  mutate(atlas = "self-referential") %>%
  bind_rows(ggeffects::ggpredict(mod_h4bm, c("cond", "topic", "pID"), type = "random") %>%
              data.frame() %>%
              mutate(atlas = "mentalizing")) %>%
  mutate(x = factor(x, levels = c("self", "control", "other")),
         atlas = factor(atlas, levels = c("self-referential", "mentalizing")))

predicted_h4m %>%
  ggplot(aes(x = x, y = predicted, color = group)) +
  stat_summary(data = predicted_sub_h4m, aes(group = interaction(group, facet)), fun = "mean", geom = "line", size = .1) +
  stat_summary(aes(group = group), fun = "mean", geom = "line", size = 1, position = position_dodge(.1)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high, group = group),
                  size = .75, position = position_dodge(.1)) +
  facet_grid(~atlas) +
  scale_color_manual(name = "", values = palette_topic) +
  labs(x = "", y = "ROI activity (SD)\n") +
  plot_aes +
  theme(legend.position = c(.18, .95))

H5

Are the effect of the experimental manipulations on sharing intentions moderated by article topic?

There is a main effect of topic, such that health articles have higher sharing intentions than climate articles.

These data are not consistent with moderation by topic.

mod_h5m = lmer(value ~ cond * topic + (1 + cond + topic | pID),
              data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

plot

predicted_h5m = ggeffects::ggpredict(mod_h5m, c("cond", "topic")) %>%
  data.frame() %>%
  mutate(x = factor(x, levels = c("self", "control", "other")))

predicted_sub_h5m = ggeffects::ggpredict(mod_h5m, terms = c("cond", "topic", "pID"), type = "random") %>%
  data.frame() %>%
  mutate(x = factor(x, levels = c("self", "control", "other")))
  
predicted_h5m %>%
  ggplot(aes(x = x, y = predicted, color = group)) +
  stat_summary(data = predicted_sub_h5m, aes(group = interaction(group, facet)), fun = "mean", geom = "line", size = .1) +
  stat_summary(aes(group = group), fun = "mean", geom = "line", size = 1) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high, group = group),
                  size = .75) +
  scale_color_manual(name = "", values = palette_topic) +
  labs(x = "", y = "predicted sharing intention\n") +
  plot_aes +
  theme(legend.position = c(.85, .15))

model table

table_h5m = table_model(mod_h5m)

table_h5m %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.53 [2.43, 2.62] 110.42 52.81 < .001
other -0.06 [-0.15, 0.03] 241.50 -1.30 .196
self -0.04 [-0.13, 0.04] 301.81 -0.97 .331
topic (health) 0.23 [0.10, 0.37] 146.39 3.40 < .001
other x topic (health) 0.05 [-0.07, 0.16] 5683.42 0.83 .404
self x topic (health) -0.01 [-0.12, 0.11] 5682.25 -0.09 .931

summary

summary(mod_h5m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: value ~ cond * topic + (1 + cond + topic | pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 16292.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.65946 -0.72242  0.04713  0.74522  2.59120 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr             
##  pID      (Intercept) 0.12180  0.3490                    
##           condother   0.03257  0.1805   -0.13            
##           condself    0.02347  0.1532   -0.25  0.90      
##           topichealth 0.25714  0.5071   -0.32 -0.20 -0.50
##  Residual             0.85254  0.9233                    
## Number of obs: 5935, groups:  pID, 85
## 
## Fixed effects:
##                          Estimate  Std. Error          df t value
## (Intercept)              2.529544    0.047899  110.416538  52.810
## condother               -0.059513    0.045935  241.496147  -1.296
## condself                -0.043601    0.044770  301.814917  -0.974
## topichealth              0.234436    0.068927  146.387458   3.401
## condother:topichealth    0.049004    0.058720 5683.423862   0.835
## condself:topichealth    -0.005068    0.058748 5682.245413  -0.086
##                                   Pr(>|t|)    
## (Intercept)           < 0.0000000000000002 ***
## condother                         0.196356    
## condself                          0.330885    
## topichealth                       0.000865 ***
## condother:topichealth             0.404019    
## condself:topichealth              0.931259    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cndthr cndslf tpchlt cndth:
## condother   -0.435                            
## condself    -0.474  0.562                     
## topichealth -0.465  0.205  0.129              
## cndthr:tpch  0.306 -0.640 -0.327 -0.426       
## cndslf:tpch  0.306 -0.319 -0.657 -0.426  0.500

H6

Are the relationships between ROI activity positively and sharing intentions moderated by article topic?

self-referential ROI

There is a main effect of topic, such that health articles have higher sharing intentions than climate articles.

These data are not consistent with moderation by topic.

mod_h6am = lmer(value ~ self_referential * topic + (1 + self_referential + topic | pID),
              data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h6am = table_model(mod_h6am)

table_h6am %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.49 [2.41, 2.57] 84.51 60.38 < .001
self-referential 0.06 [0.03, 0.10] 2284.13 3.72 < .001
topic (health) 0.24 [0.12, 0.36] 84.76 3.99 < .001
self-referential x topic (health) 0.02 [-0.03, 0.07] 5554.69 0.83 .405

summary

summary(mod_h6am)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: value ~ self_referential * topic + (1 + self_referential + topic |  
##     pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 16264.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.74920 -0.72685  0.05018  0.74273  2.71245 
## 
## Random effects:
##  Groups   Name             Variance Std.Dev. Corr       
##  pID      (Intercept)      0.119488 0.34567             
##           self_referential 0.000563 0.02373  -1.00      
##           topichealth      0.254015 0.50400  -0.43  0.35
##  Residual                  0.853769 0.92400             
## Number of obs: 5935, groups:  pID, 85
## 
## Fixed effects:
##                                Estimate Std. Error         df t value
## (Intercept)                     2.48872    0.04121   84.51383  60.385
## self_referential                0.06373    0.01713 2284.12966   3.720
## topichealth                     0.23877    0.05983   84.75844   3.991
## self_referential:topichealth    0.01993    0.02394 5554.68973   0.833
##                                          Pr(>|t|)    
## (Intercept)                  < 0.0000000000000002 ***
## self_referential                         0.000204 ***
## topichealth                              0.000139 ***
## self_referential:topichealth             0.405036    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) slf_rf tpchlt
## self_rfrntl -0.178              
## topichealth -0.473  0.073       
## slf_rfrntl:  0.030 -0.696 -0.058

mentalizing ROI

There is a main effect of topic, such that health articles have higher sharing intentions than climate articles.

These data are not consistent with moderation by topic.

mod_h6bm = lmer(value ~ mentalizing * topic + (1 + topic | pID),
              data = merged_wide,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h6bm = table_model(mod_h6bm)

table_h6bm %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 2.47 [2.39, 2.55] 87.14 59.32 < .001
mentalizing 0.07 [0.04, 0.10] 5544.20 4.18 < .001
topic (health) 0.25 [0.13, 0.37] 88.06 4.07 < .001
mentalizing x topic (health) -0.01 [-0.06, 0.04] 5609.50 -0.37 .713

summary

summary(mod_h6bm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: value ~ mentalizing * topic + (1 + topic | pID)
##    Data: merged_wide
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 16275.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.75119 -0.72361  0.05105  0.75218  2.66333 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  pID      (Intercept) 0.1204   0.3469        
##           topichealth 0.2551   0.5050   -0.43
##  Residual             0.8554   0.9249        
## Number of obs: 5935, groups:  pID, 85
## 
## Fixed effects:
##                            Estimate  Std. Error          df t value
## (Intercept)                2.472130    0.041672   87.144707  59.323
## mentalizing                0.070210    0.016777 5544.196558   4.185
## topichealth                0.246283    0.060515   88.059778   4.070
## mentalizing:topichealth   -0.008819    0.023953 5609.499240  -0.368
##                                     Pr(>|t|)    
## (Intercept)             < 0.0000000000000002 ***
## mentalizing                         0.000029 ***
## topichealth                         0.000102 ***
## mentalizing:topichealth             0.712762    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mntlzn tpchlt
## mentalizing -0.133              
## topichealth -0.482  0.091       
## mntlzng:tpc  0.092 -0.695 -0.149

combined plot

vals = seq(-4.5,4.5,.1)

predicted_h6m = ggeffects::ggpredict(mod_h6am, c("self_referential [vals]", "topic")) %>%
  data.frame() %>%
  mutate(atlas = "self-referential") %>%
  bind_rows(ggeffects::ggpredict(mod_h6bm, c("mentalizing [vals]", "topic")) %>%
              data.frame() %>%
              mutate(atlas = "mentalizing")) %>%
  mutate(atlas = factor(atlas, levels = c("self-referential", "mentalizing")))

predicted_sub_h6m = ggeffects::ggpredict(mod_h6am, terms = c("self_referential [vals]", "topic", "pID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "self-referential") %>%
  bind_rows(ggeffects::ggpredict(mod_h6bm, c("mentalizing [vals]", "topic", "pID"), type = "random") %>%
              data.frame() %>%
              mutate(roi = "mentalizing")) %>%
  mutate(roi = factor(roi, levels = c("self-referential", "mentalizing")))

predicted_h6m %>%
  ggplot(aes(x = x, y = predicted, color = group, fill = group)) +
  stat_smooth(data = predicted_sub_h6m, aes(group = interaction(group, facet)),
              geom ='line', method = "lm", alpha = .1, size = 1, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .3, color = NA) +
  geom_line(size = 2) +
  facet_grid(~atlas) +
  scale_color_manual(name = "", values = palette_topic) +
  scale_fill_manual(name = "", values = palette_topic) +
  labs(y = "predicted sharing intention\n", x = "\nROI activity (SD)") +
  plot_aes +
  theme(legend.key.width=unit(2,"cm"))

combined table

table_h1am %>% mutate(DV = "H1a: Self-relevance") %>%
  bind_rows(table_h1bm %>% mutate(DV = "H1b: Social relevance")) %>%
  bind_rows(table_h2am %>% mutate(DV = "H2a: Self-relevance")) %>%
  bind_rows(table_h2bm %>% mutate(DV = "H2b: Social relevance")) %>%
  bind_rows(table_h3m %>% mutate(DV = "H3a-b: Sharing intention")) %>%
  bind_rows(table_h4am %>% mutate(DV = "H4a: Self-referential ROI")) %>%
  bind_rows(table_h4bm %>% mutate(DV = "H4b: Mentalizing ROI")) %>%
  bind_rows(table_h5m %>% mutate(DV = "H5: Sharing intention")) %>%
  bind_rows(table_h6am %>% mutate(DV = "H6a: Sharing intention")) %>%
  bind_rows(table_h6bm %>% mutate(DV = "H6b: Sharing intention")) %>%
  select(DV, everything()) %>%
  kable() %>%
  kable_styling()
DV term b [95% CI] df t p
H1a: Self-relevance intercept 2.49 [2.40, 2.59] 83.88 51.73 < .001
H1a: Self-relevance self-referential 0.03 [-0.01, 0.06] 5807.44 1.56 .120
H1a: Self-relevance topic (health) 0.13 [0.01, 0.26] 84.45 2.14 .035
H1a: Self-relevance self-referential x topic (health) 0.04 [-0.01, 0.08] 5759.32 1.49 .137
H1b: Social relevance intercept 2.51 [2.43, 2.60] 97.93 58.12 < .001
H1b: Social relevance mentalizing 0.03 [0.00, 0.07] 226.60 2.16 .032
H1b: Social relevance topic (health) 0.29 [0.24, 0.34] 5924.99 12.18 < .001
H1b: Social relevance mentalizing x topic (health) 0.01 [-0.03, 0.05] 5903.41 0.66 .512
H2a: Self-relevance intercept 2.44 [2.34, 2.55] 128.84 45.55 < .001
H2a: Self-relevance other 0.04 [-0.04, 0.12] 5840.55 1.03 .304
H2a: Self-relevance self 0.12 [0.04, 0.20] 5840.64 2.84 .004
H2a: Self-relevance topic (health) 0.22 [0.08, 0.36] 139.35 3.18 .002
H2a: Self-relevance other x topic (health) -0.07 [-0.19, 0.04] 5840.57 -1.28 .202
H2a: Self-relevance self x topic (health) -0.17 [-0.28, -0.06] 5840.67 -2.92 .004
H2b: Social relevance intercept 2.49 [2.39, 2.60] 122.27 47.37 < .001
H2b: Social relevance other 0.03 [-0.04, 0.10] 5840.82 0.78 .437
H2b: Social relevance self 0.07 [-0.00, 0.15] 5840.90 1.95 .051
H2b: Social relevance topic (health) 0.31 [0.19, 0.43] 151.89 5.04 < .001
H2b: Social relevance other x topic (health) 0.03 [-0.07, 0.14] 5840.82 0.58 .561
H2b: Social relevance self x topic (health) -0.06 [-0.16, 0.05] 5840.94 -1.08 .282
H3a-b: Sharing intention intercept 2.56 [2.50, 2.61] 113.93 86.74 < .001
H3a-b: Sharing intention self-relevance 0.30 [0.24, 0.35] 275.19 10.76 < .001
H3a-b: Sharing intention topic (health) 0.14 [0.10, 0.19] 5544.95 6.43 < .001
H3a-b: Sharing intention social relevance 0.20 [0.14, 0.26] 180.07 6.53 < .001
H3a-b: Sharing intention self-relevance x topic (health) 0.03 [-0.03, 0.09] 5555.76 1.04 .297
H3a-b: Sharing intention topic (health) x social relevance 0.06 [0.00, 0.12] 5628.34 1.99 .047
H4a: Self-referential ROI intercept 0.01 [-0.10, 0.13] 101.71 0.24 .812
H4a: Self-referential ROI other 0.09 [-0.01, 0.19] 222.78 1.76 .080
H4a: Self-referential ROI self 0.13 [0.02, 0.23] 193.75 2.38 .018
H4a: Self-referential ROI topic (health) 0.14 [0.05, 0.23] 528.65 2.99 .003
H4a: Self-referential ROI other x topic (health) -0.01 [-0.13, 0.12] 5677.79 -0.09 .928
H4a: Self-referential ROI self x topic (health) -0.08 [-0.20, 0.05] 5678.30 -1.22 .221
H4b: Mentalizing ROI intercept 0.27 [0.16, 0.39] 102.44 4.63 < .001
H4b: Mentalizing ROI other 0.06 [-0.04, 0.15] 233.15 1.16 .249
H4b: Mentalizing ROI self 0.10 [-0.00, 0.20] 198.21 1.91 .057
H4b: Mentalizing ROI topic (health) 0.11 [0.03, 0.20] 606.69 2.52 .012
H4b: Mentalizing ROI other x topic (health) 0.01 [-0.12, 0.13] 5677.90 0.09 .927
H4b: Mentalizing ROI self x topic (health) -0.06 [-0.18, 0.07] 5678.40 -0.89 .373
H5: Sharing intention intercept 2.53 [2.43, 2.62] 110.42 52.81 < .001
H5: Sharing intention other -0.06 [-0.15, 0.03] 241.50 -1.30 .196
H5: Sharing intention self -0.04 [-0.13, 0.04] 301.81 -0.97 .331
H5: Sharing intention topic (health) 0.23 [0.10, 0.37] 146.39 3.40 < .001
H5: Sharing intention other x topic (health) 0.05 [-0.07, 0.16] 5683.42 0.83 .404
H5: Sharing intention self x topic (health) -0.01 [-0.12, 0.11] 5682.25 -0.09 .931
H6a: Sharing intention intercept 2.49 [2.41, 2.57] 84.51 60.38 < .001
H6a: Sharing intention self-referential 0.06 [0.03, 0.10] 2284.13 3.72 < .001
H6a: Sharing intention topic (health) 0.24 [0.12, 0.36] 84.76 3.99 < .001
H6a: Sharing intention self-referential x topic (health) 0.02 [-0.03, 0.07] 5554.69 0.83 .405
H6b: Sharing intention intercept 2.47 [2.39, 2.55] 87.14 59.32 < .001
H6b: Sharing intention mentalizing 0.07 [0.04, 0.10] 5544.20 4.18 < .001
H6b: Sharing intention topic (health) 0.25 [0.13, 0.37] 88.06 4.07 < .001
H6b: Sharing intention mentalizing x topic (health) -0.01 [-0.06, 0.04] 5609.50 -0.37 .713

cite packages

report::cite_packages()
##   - Angelo Canty and Brian Ripley (2021). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-28.
##   - Douglas Bates, Martin Maechler and Mikael Jagan (2023). Matrix: Sparse and Dense Matrix Classes and Methods. R package version 1.5-4. https://CRAN.R-project.org/package=Matrix
##   - Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
##   - H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
##   - Hadley Wickham (2021). forcats: Tools for Working with Categorical Variables (Factors). R package version 0.5.1. https://CRAN.R-project.org/package=forcats
##   - Hadley Wickham (2022). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.0. https://CRAN.R-project.org/package=stringr
##   - Hadley Wickham and Maximilian Girlich (2022). tidyr: Tidy Messy Data. R package version 1.2.0. https://CRAN.R-project.org/package=tidyr
##   - Hadley Wickham, Jennifer Bryan and Malcolm Barrett (2021). usethis: Automate Package and Project Setup. R package version 2.1.5. https://CRAN.R-project.org/package=usethis
##   - Hadley Wickham, Jim Hester and Jennifer Bryan (2022). readr: Read Rectangular Text Data. R package version 2.1.2. https://CRAN.R-project.org/package=readr
##   - Hadley Wickham, Jim Hester, Winston Chang and Jennifer Bryan (2021). devtools: Tools to Make Developing R Packages Easier. R package version 2.4.3. https://CRAN.R-project.org/package=devtools
##   - Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2022). dplyr: A Grammar of Data Manipulation. R package version 1.0.9. https://CRAN.R-project.org/package=dplyr
##   - Hao Zhu (2021). kableExtra: Construct Complex Table with 'kable' and Pipe Syntax. R package version 1.3.4. https://CRAN.R-project.org/package=kableExtra
##   - Jim Hester, Hadley Wickham and Gábor Csárdi (2021). fs: Cross-Platform File System Operations Based on 'libuv'. R package version 1.5.2. https://CRAN.R-project.org/package=fs
##   - Kirill Müller and Hadley Wickham (2022). tibble: Simple Data Frames. R package version 3.1.8. https://CRAN.R-project.org/package=tibble
##   - Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest Package:Tests in Linear Mixed Effects Models." _Journal of StatisticalSoftware_, *82*(13), 1-26. doi: 10.18637/jss.v082.i13 (URL:https://doi.org/10.18637/jss.v082.i13).
##   - Lionel Henry and Hadley Wickham (2020). purrr: Functional Programming Tools. R package version 0.3.4. https://CRAN.R-project.org/package=purrr
##   - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects fromRegression Models." _Journal of Open Source Software_, *3*(26), 772.doi: 10.21105/joss.00772 (URL: https://doi.org/10.21105/joss.00772).
##   - R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
##   - Rinker, T. W. & Kurkiewicz, D. (2017). pacman: Package Management for R. version 0.5.0. Buffalo, New York. http://github.com/trinker/pacman
##   - Wickham et al., (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686
##   - Yihui Xie (2021). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.37.